%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SVM Decoding
%
% GE_Step1b_plot_individual.m
%
% Modified from Dimitrios Pantazis's code 
%  by SA 2018-03-30, 2018-04-11
%
% SA 2019-08-09: modified from GE_Step1_classifierapply.m
%                to make nicely formatted plots of 
%                the classifier accuracy for all condition pairs
%                for each ROI and each subject
%
%  Input: *_all.mat file from "results" subdirectory for each subject
%
%  Output: *_accuracy.fig and *accuracy.png files 
%          in "figures" subdirectory for each subject

%initialize

clear all;
close all;

svm_dir = '/autofs/space/clive_001/users/adriana/GE_SVM';
scripts_dir = sprintf('%s/scripts', svm_dir);
%functions_dir = sprintf('%s/scripts/Functions', svm_dir);

%cd(svm_dir);

addpath(genpath(scripts_dir));
%addpath(genpath(functions_dir));

%SubjectNames = {'GE_05', 'GE_06', 'GE_07', 'GE_08', 'GE_09', 'GE_10', 'GE_11', 'GE_12', 'GE_13', 'GE_14', 'GE_15', 'GE_16'};
SubjectNames = {'GE_12', 'GE_13'}; %{'GE_05', 'GE_06', 'GE_07', 'GE_08', 'GE_09', 'GE_10', 'GE_11', 'GE_12', 'GE_13', 'GE_14', 'GE_15', 'GE_16', 'GE_17', 'GE_18', 'GE_19', 'GE_20', 'GE_21', 'GE_22', 'GE_23', 'GE_24'};
n_subjects = length(SubjectNames);

%analysis_tag = '_neighbors_vs_seeds_16Slices_2Rois_SMG_pMTG'
analysis_tag = ''

% Specify the Regions of Interest ("ROIs", also called "labels":
%roiset_array = {'pMTG'};
%roiset_array = {'LIFG', 'SMG'};
%roiset_array = {'FrontalPole', 'CentGyri'};
roiset_array = {'L_ITG_1-lh','L_ITG_2-lh', 'L_MTG_1-lh', 'L_MTG_2-lh', 'L_ParaHip_1-lh', 'L_ParsTri_1-lh', 'L_postCG_1-lh', 'L_SMG_1-lh', 'L_STG_1-lh', 'R_ITG_1-rh', 'R_ITG_2-rh', 'R_MTG_1-rh', 'R_MTG_2-rh', 'R_MTG_3-rh', 'R_MTG_4-rh', 'R_postCG_2-rh', 'R_postCG_4-rh', 'R_preCG_1-rh', 'R_preCG_2-rh', 'R_STG_1-rh', 'R_STG_2-rh', 'R_STG_3-rh'};
n_roisets = length(roiset_array);

% Specify the condition pairs:
condition_tag = 'nvs'; % 'nvs'|'lvn'|'pos'

switch condition_tag
	case 'nvs'
		%Train all neighbors; test seeds
		categories = {'neighbors'};
		train_conditionA_array = {'11','11','11','11','11','22','22','22','22','33','33','33','44','44','55'};
		train_conditionB_array = {'22','33','44','55','66','33','44','55','66','44','55','66','55','66','66'};
		test_conditionA_array  = { '1', '1', '1', '1', '1', '2', '2', '2', '2', '3', '3', '3', '4', '4', '5'};
		test_conditionB_array  = { '2', '3', '4', '5', '6', '3', '4', '5', '6', '4', '5', '6', '5', '6', '6'};
	case 'lvn'
		%Separately train lexical and nonword neighbors; test seeds
		categories = {'nonwords'; 'words'};
		train_conditionA_array = {'100','100','100','100','100','200','200','200','200','300','300','300','400','400','500';...
									'10','10','10','10','10','20','20','20','20','30','30','30','40','40','50'};
		train_conditionB_array = {'200','300','400','500','600','300','400','500','600','400','500','600','500','600','600';...
									'20','30','40','50','60','30','40','50','60','40','50','60','50','60','60'};
		test_conditionA_array = {'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5';...
									'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5'};
		test_conditionB_array = {'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6';...
									'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6'};
	case 'pos'
		%Separately train based on differential position; test seeds
		categories = {'initial'; 'vowel'; 'final'};
		train_conditionA_array = {'1001', '1001', '1001', '1001', '1001', '2001', '2001', '2001', '2001', '3001', '3001', '3001', '4001', '4001', '5001';...
									'1002', '1002', '1002', '1002', '1002', '2002', '2002', '2002', '2002', '3002', '3002', '3002', '4002', '4002', '5002';...
									'1003', '1003', '1003', '1003', '1003', '2003', '2003', '2003', '2003', '3003', '3003', '3003', '4003', '4003', '5003'};
		train_conditionB_array = {'2001', '3001', '4001', '5001', '6001', '3001', '4001', '5001', '6001', '4001', '5001', '6001', '5001', '6001', '6001';...
									'2002', '3002', '4002', '5002', '6002', '3002', '4002', '5002', '6002', '4002', '5002', '6002', '5002', '6002', '6002';...
									'2003', '3003', '4003', '5003', '6003', '3003', '4003', '5003', '6003', '4003', '5003', '6003', '5003', '6003', '6003'};
		test_conditionA_array = {'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5';...
									'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5';...
									'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5'};
		test_conditionB_array = {'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6';...
									'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6';...
									'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6'};
end

n_condpairs = min(length(train_conditionA_array), length(train_conditionB_array));

% These are for plotting with subplots:
n_hubwords = 6;
n_cols = n_hubwords - 1;
n_rows = n_hubwords - 1;

% Select the time window for linear fit (in milliseconds):
%fit_tmin = -100;
%fit_tmax = 1000;
fit_tmin = 0;
fit_tmax = 800;
%fit_tmin = 100;
%fit_tmax = 900;

% plot_only = 0;   % compute and visualize
plot_only = 1;  % do not re-compute, but only visualize the results

%my_file_tag = '_seppo';
my_file_tag = '';

%-----------------------------------------------------------------
for i_category=1:length(categories)

for i_subject = 1:n_subjects %for all subjects
    %i_subject = 1;
    subject = char(SubjectNames(i_subject));
    %subject_dir=char(SubjectDirectories(i_subject));
    subject_dir = char(strcat(subject, analysis_tag));
    
    % Create the "<subject>/figures" subdirectory in case it doesn't exist yet
    figure_dir = sprintf('%s/%s/figures', svm_dir, subject_dir);
    [status, msg, msgID] = mkdir(figure_dir);

    for i_roiset = 1:n_roisets
        roiset = char(roiset_array(i_roiset));
        
        figure('Position',[0, 0, 1650, 1200]);
        %figure(i_subject*100+i_roiset); % separate figure for each ROI (and each subject)
        
        for i_condpair = 1:n_condpairs
        
            train_conditionA= char(train_conditionA_array(i_category, i_condpair));
            train_conditionB= char(train_conditionB_array(i_category, i_condpair));
            test_conditionA= char(test_conditionA_array(i_category, i_condpair));
            test_conditionB= char(test_conditionB_array(i_category, i_condpair));

            %plot previously computed SVM results:

            resultsfile = sprintf('%s/%s/results/%s_%s_train%sand%s_test%svs%s_Accuracy.mat', svm_dir, subject_dir, subject, roiset,train_conditionA,train_conditionB,test_conditionA,test_conditionB);
            load(resultsfile)
            if exist('x')
                accuracy = x;
                Time = y;
            end
            % Time window for linear fit:
            % sfreq = 1000.;
            %sfreq = 1./(Time(2)-Time(1));  % sfreq should be taken into
            %account below when converting fit_tmin and fit_tmax to fit_i1
            %and fit_i2; however, since in this study it is always 1 (kHz),
            %we ignore it in the formula
            fit_i1 = fit_tmin - round(Time(1));
            fit_i2 = fit_tmax - round(Time(1));

            fitTime = Time(fit_i1:fit_i2);
            %fitAccuracy = accuracy(i_subject,fit_i1:fit_i2);
            fitAccuracy = accuracy(1,fit_i1:fit_i2);
            fitobject = fit(fitTime',mean(fitAccuracy,1)','poly1');
            %fitobject = polyfit(fitTime',mean(fitAccuracy,1)',1);
            %fitfile = sprintf('%s/%s/results/%s_%s_train%sand%s_test%svs%sfitline.mat', svm_dir, subject_dir, subject, roiset, train_conditionA, train_conditionB, test_conditionA, test_conditionB);
            %if (~plot_only)
            %    save(fitfile,'fitobject');
            %end
              
            iA = round(str2double(test_conditionA));
            iB = round(str2double(test_conditionB));
            subplot(n_rows, n_cols, (iA-1)*n_cols +  iB-1);
            plot(fitobject, fitTime, mean(fitAccuracy,1));
            %fitobject = fit(Time',mean(accuracy,1)','poly1');
            %plot(fitobject, Time, mean(accuracy,1));
            set(gca,'fontsize',8);
            ylabel('Accuracy (%)','fontsize',10)
            xlabel('Time (msec)','fontsize',10)
            yline(50);

            %titletext = sprintf('%s %s: train %s/%s, test %s/%s', subject_dir, roiset, train_conditionA, train_conditionB,test_conditionA, test_conditionB);
            maxlength = 33; % cutoff for too long title information 
            if length(subject_dir) > maxlength
                plot_title = subject_dir(1:maxlength);
            else
                plot_title = subject_dir;
            end
            titletext = sprintf('%s %s: train %s/%s, test %s/%s', plot_title, roiset, train_conditionA, train_conditionB,test_conditionA, test_conditionB);
            title(titletext, 'fontsize', 10,'Interpreter','none');
            axis([-200 1100 0 100]);
            legend('off');

            %results_fig_file = sprintf('%s/%s/results/%s_%s_%svs%s_AverageAccuracy.png', svm_dir, subject, subject, roiset, condA, condB)
            %print(gcf,[results_fig_file],'-dpng','-r300');

            % Store decoding accuracy for all condition pairs for this subject:
            accuracy_allpairs(i_condpair,:) = accuracy; 

        end % i_condpair
    
        % Mean accuracy over all condition pairs:

        accuracy_ave = mean(accuracy_allpairs,1);
        accuracy_std = std(accuracy_allpairs,1);
        accuracy_sem = accuracy_std/sqrt(n_condpairs);

        % Plot the mean accuracy in the lower left corner of the figure:
        
        subplot(n_rows, n_cols, (n_rows-1)*n_cols + 1);
        plot(Time, accuracy_ave);
        shadedErrorBar(Time,accuracy_ave,accuracy_sem,'lineProps','-r','patchSaturation',0.16);
        set(gca,'fontsize',8);
        ylabel('Accuracy (%)','fontsize',10)
        xlabel('Time (msec)','fontsize',10)
        yline(50);

        plot_title = subject_dir;
        maxlength = 33; % cutoff for too long title information 
        if length(subject_dir) > maxlength
            plot_title = subject_dir(1:maxlength);
        end
        
        titletext = sprintf('%s %s: all cond. pairs, %s', plot_title, roiset, categories{i_category});
        title(titletext, 'fontsize', 10,'Interpreter','none');
        axis([-200 1100 0 100]);
        legend('off');
        
        % Save the figure:
        %save the average accuracy across pairwise conditions as json
        category_dir = sprintf('%s/%s/results/%s', svm_dir, subject, categories{i_category});
        mkdir(category_dir);
        %%%temporary workaround for poor filenaming
        %format_roiset = [strtok(roiset, '1'), '_1']
        %roiset = format_roiset;
        %%%
		format_roiset = strtok(roiset, '-');
        ave_fname = sprintf('%s/%s_%s_%s_aveAccuracy.mat', category_dir, subject, format_roiset, categories{i_category});
		category = categories{i_category};
		save(ave_fname, 'subject', 'format_roiset', 'category', 'accuracy_ave', 'Time');
        %json_struct = struct("Subject",subject, "ROI",roiset, "Category", categories{i_category}, "Average_Accuracy",accuracy_ave,"Time",Time);
        %encoded = jsonencode(json_struct);
        %fid = fopen(ave_fname, 'w');
        %fprintf(fid, encoded);
        %fclose(fid);

        %fig_file = sprintf('%s/%s/results/%s_%d-%dms_AverageAccuracy%s.fig', svm_dir, subject_dir, subject, fit_tmin, fit_tmax, my_file_tag)
        fig_file = sprintf('%s/%s_%s_%d-%dms_accuracy_%s.fig', figure_dir, subject, roiset, fit_tmin, fit_tmax, categories{i_category})
        %savefig(fig_file);

        %png_file = sprintf('%s/%s/results/%s_%d-%dms_AverageAccuracy%s.png', svm_dir, subject_dir, subject, fit_tmin, fit_tmax, my_file_tag)
        png_file = sprintf('%s/%s_%s_%d-%dms_accuracy_%s.png', figure_dir, subject, roiset, fit_tmin, fit_tmax, categories{i_category})
        %print(gcf,[png_file],'-dpng','-r300');

    end % i_roiset

end % i_subject

end % i_category




